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ABSTRACT 

We extend our models of the vertical structure and emergent radiation field of accretion disks around super- 
massive black holes described in previous papers of this series. Our models now include both a self-consistent 
treatment of Compton scattering and the effects of continuum opacities of the most important metal species (C, 
N, O, Ne, Mg, Si, S, Ar, Ca, Fe, Ni). With these new effects incorporated, we compute the predicted spectrum 
from black holes accreting at nearly the Eddington luminosity (L/LEdd ~ 0.3) and central masses of 10'', 10^, and 
1O^M0. We also consider two values of the Shakura-Sunyaev a parameter, 0.1 and 0.01, but in contrast to our 
previous papers, we consider a kinematic viscosity which is independent of depth. 

Although it has little effect when M > lO^M©, Comptonization grows in importance as the central mass de- 
creases and the central temperature rises. It generally produces an increase in temperature with height in the 
uppermost layers of hot atmospheres. Compared to models with coherent electron scattering, Comptonized mod- 
els have enhanced extreme ultraviolet/soft X-ray emission, but they also have a more sharply declining spectrum 
at very high frequencies. Comptonization also smears the hydrogen and the He II Lyman edges. The effects of 
metals on the overall spectral energy distribution are smaller than the effects of Comptonization for these param- 
eters. Compared to pure hydrogen-helium models, models with metal continuum opacities have reduced flux in 
the high frequency tail, except at the highest frequencies, where the flux is very low. Metal photoionization edges 
are not present in the overall disk-integrated model spectra. 

The viscosity parameter a has a more dramatic effect on the emergent spectrum than do metal continuum opac- 
ities. As a increases (and therefore the disk column density decreases), the flux at both the high and low frequency 
extremes of the spectrum increases, while the flux near the peak decreases. Multitemperature blackbodies are a 
very poor approximation to accretion disk spectra in the soft X-ray region, and such crude modeling may greatly 
overestimate the accretion luminosity required to explain observed soft X-ray excesses in active galactic nuclei. 
In addition to our new grid of models, we also present a simple analytic prescription for the vertical temperature 
structure of the disk in the presence of Comptonization, and show under what conditions a hot outer layer (a 
corona) is formed. 

Subject headings: accretion, accretion disks — galaxies; active — galaxies:nuclei — radiative transfer 



1. INTRODUCTION 

In previous papers of this series (Hubeny & Hubeny 1997, 
1998; Hubeny et al. 2000; hereafter referred to as Papers I, II, 
and III, respectively), we computed a series of models of ac- 
cretion disks around supermassive black holes appropriate for 
quasars and active galactic nuclei (AGN). The models included 
a set of 1-D vertical structure calculations for the individual 
annuli of the disk. We solved self-consistently the equations 
of vertical hydrostatic equilibrium, energy balance, radiative 
transfer, and, since we did not assume local thermodynamic 
equilibrium (LTE), the set of statistical equilibrium equations 
for selected energy levels of various species of atoms and ions. 
In these previous models, we considered a simple H-He chem- 
ical composition, and we assumed that electron scattering was 
coherent (i.e., Thomson scattering). 

'present address: California Institute of Technology, Pasadena, CA 91125 



In the present paper, we refine the previous modeling proce- 
dure by relaxing these two limiting approximations. We now in- 
corporate metals, and we include the effects of incoherent elec- 
tron scattering, i.e. Compton scattering. Both of these effects 
are expected to modify the emergent spectrum of individual 
hot annuli which emit predominantly in the extreme ultraviolet 
and soft X-rays. Excess emission above a power law extending 
down from the hard X-rays is widely observed in many classes 
of AGN, from Seyferts to quasars (e.g. Turner & Pounds 1989). 
The origin of this emission remains unclear, but it has been pro- 
posed to be thermal, Comptonized emission from hot portions 
of an accretion disk (e.g. Czerny & Elvis 1987; Ross, Fabian & 
Mineshige 1992). In addition, fluorescent line emission and 
Compton reflection humps are observed in AGN (Nandra & 
Pounds 1994) and appear to be at least partly the result of exter- 
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nal illumination of the disk by hard X-rays from a hot, coronal 
region somewhere in the vicinity. The effects of Compton scat- 
tering and metal opacities within the disk itself are crucial for 
understanding this phenomenology. Moreover, these effects are 
still more pronounced in the even higher temperature accretion 
disks thought to exist around stellar mass black hole candidates 
in X-ray binaries. 

This paper is organized as follows. In section 2 we describe 
in detail our basic physical assumptions and numerical methods 
in treating the effects of Comptonization and metal opacities. 
We then describe in section 3 an approximate analytic calcu- 
lation of the vertical temperature structure of an atmosphere 
which is subject to Comptonization. In section 4 we present 
our new numerical models of hot accretion disks, including a 
detailed comparison with our previous work and that of other 
authors. In section 5 we discuss these results and then we sum- 
marize our conclusions in section 6. 

2. METHOD 

The basic equations for the vertical structure were discussed 
in detail in Paper II, to which the reader is referred for more 
information. The main difference here is an inclusion of Comp- 
ton scattering. Since our treatment differs somewhat from more 
commonly used approaches (e.g., Ross, Fabian, & Mineshige 
1992; Shimura & Takahai-a 1993 in the context of AGN disks; 
or Madej 1989, 1991, 1998 in the context of model atmospheres 
of neutron stars and white dwarfs), we will describe it here in 
more detail. 

2.1. Treatment of Compton Scattering 

A general treatment of Compton scattering of (unpolar- 
ized) radiation is described in Pomraning (1973). The ef- 
fects of Compton scattering are described through an appro- 
priate frequency- and angle-dependent redistribution function 
that incorporates Klein-Nishina corrections and the relativistic 
Maxwellian velocity distribution of the electrons. In the context 
of stellar atmosphere modeling, this approach was used e.g. by 
Madej (1989; 1991, 1998; Madej & Rozanska 2000) for com- 
puting model atmospheres of neutron stars, white dwarfs, and 
X-ray irradiated B stars. 

We are interested here in accurately treating the extreme ul- 
traviolet/soft X-ray emission of annuli with temperatures much 
less than the electron rest mass energy. In this non-relativistic 
regime, Compton scattering can be incorporated in the radia- 
tive transfer equation through a Kompaneets-like term, i.e. by 
replacing the exact angle-frequency redistribution integral by 
a Fokker-Planck form. We present details of how this can be 
implemented numerically in Appendix A. For all the calcula- 
tions presented in this paper, we use a Compton source func- 
tion which is averaged over all angles, analogous to the stan- 
dard isotropic Kompaneets equation (e.g. Rybicki & Lightman 
1979). We stress, however, that we do not assume that the 
radiation field is isotropic. Once an atmosphere structure has 
converged, our treatment could be improved by performing a 
formal solution of the radiative transfer equation with the full 
angle-dependent Compton source function (see Appendix A). 

We adopt the notation commonly used in stellar atmosphere 
modeling, i.e. we express a photon's energy through its fre- 
quency, V, and the photon occupation number through the spe- 
cific intensity of radiation, Also, we often formally express 
a dependence on frequency through a subscript v. In the plane- 
parallel atmospheres we are considering, the radiative transfer 



equation may be written as 



(1) 



Here ^ is the direction cosine between the photon propagation 
direction and the vertical. The monochromatic optical depth, 
r^, is defined by 



(2) 



where is the total thermal absorption coefficient at frequency 
V (including contributions from all bound-free, free-free, and 
possibly bound-bound transitions which operate at this fre- 
quency), and the scattering coefficient is given by 



(Ty = ngCTT ( 1 — 2x) . 



(3) 



Here crj is the Thomson cross-section, and the dimensionless 
frequency x is given by 



hv 



(4) 



The — 2x term in equation ^ accounts for the Klein-Nishina 
reduction in scattering cross-section, to first order in x. 
We define the "photon destruction probability", e^, as 



(5) 



and the "scattering probability", Ajy, as 
n^cTi 1 - 



(6) 



Note that in contrast to the case of coherent scattering, X^, ^ 
1 — e^. The thermal source function is given by 



th 

5* = ^ , (7) 



where 77* is the thermal emission coefficient, which contains 
contributions from all the spontaneous inverse processes to 
those giving rise to thermal absorption at frequency v. As is 
customary, we treat the stimulated emission for these processes 
as negative absorption. 

The angle- averaged "Compton source function" is given by 



^compt ^ ( ^ _ + (x - 36)7^ + QJ'l + —-JplxiJl - J,) 

2hv^ 



where 



kT 



(8) 
(9) 



Jp is the angle-averaged intensity, and J'^ and 7" are defined as 

(10) 



rlJ 8^1 
J„ = — and J„ 



d\\\i 



d{\nvf 



A detailed derivation of equation (8) is presented in Appendix 
A, but we note that it is simply equivalent to taking the emis- 
sion terms of the isoti-opic Kompaneets equation (e.g. Rybicki 
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& Lightman 1979), and replacing the photon occupation num- 
ber with Jy/v^ . The last term on the right-hand side of equa- 
tion represents stimulated Compton scattering. Note that 
the Thomson limit of the scattering source function can be for- 
mally obtained by setting x = 8 = in equation (|]), in which 
case we obtain 5'[;°'"'" = Ju, i.e., the usual Thomson scattering 
source function for isotropic scattering. 

We write the energy balance equation in a manner very sim- 
ilar to equation (36) of Paper 11, viz. 



dF- 



= 47r 



(11) 



Here is the vertical component of the energy flux. Note that 
in the case of coherent (Thomson) scattering, the absorption 
and emission contributions exactly cancel in the right hand side 
of equation (11), so the net radiative cooling term contains only 
thermal absorption and emission terms. Compton scattering, 
however, affects the thermal balance, so we have to retain the 
total absorption coefficient Xu and the total source function, 

otot — , Cth I \ c'-ompt 

The coupled set of structural equations (|l]), (|l]), together 
with the hydrostatic equilibrium and statistical equilibrium 
equations (given explicitly in Paper 11), discretized in fre- 
quency and depth, is solved by the hybrid Complete Lineariza- 
tion/Accelerated Lambda Iteration (CL/ALl) method (Hubeny 
& Lanz 1995). The only new ingredient here, the Compton 
source function, is given by 

^Compt +B,J! + QJi+i +Ji{UiJi-i + V;/,- + W,J,-+ 1 ) , 

(12) 

where / is the frequency index, 7,- = /(t',) , and the coefficients 
A - C and f/ - W are functions of the parameters x and O. They 
also depend on the form of the difference formula for differ- 
ential operators in frequency. We use here a simple second 
order difference formula. Detailed expressions for the coeffi- 
cients are given in Appendix A. Equation (|l2|) is valid for the 
interior points of the frequency grid: i ~ 2,. . . ,NF — 1 (NF 
being the total number of frequency points). For boundary con- 
ditions in frequency space, we use the following procedure: (i) 
For the lowest frequency (/ — 1), we assume that the mean in- 
tensity follows the Planck function in the Rayleigh- Jeans limit, 
Ji, oc {2kT /c^)v^ so that 



V —11 j" — A.J 



(13) 



(ii) Analogously, for the highest frequency (/ = NF), we assume 
the Wien limit, 



2hv^ , , , 
Jv oc — T— e.x^{—hv kT) , 



(14) 



and thus 



7: = (3-z)y,, /;' = [(3-z)2-z]7,, (15) 



where 



z^hu/kT ^x/e . 



(16) 



The off-diagonal terms A,, C,, C/,, and W, therefore vanish for 
i = 1 and i = NF. 

Since the boundary condition (|l4|) is somewhat artificial, we 
have also tested another type of boundary condition, assuming 



that the mean intensity of radiation at the highest frequency is 
given by a power-law. 



in which case 



J y (3 J , J y [3 Jy 



(17) 



(18) 



We have found that the choice of boundary conditions only 
marginally influences the radiation intensity. This is because 
we adopt the highest frequency to be so high that the emergent 
flux is lower than the flux at the peak by several orders of mag- 
nitude. 

The hybrid CL/ALl scheme allows for an arbitrary parti- 
tioning of "linearized" and 'ALT' frequency points, i.e., those 
where the mean intensity of radiation is linearized, and those 
treated through the ALl scheme (and therefore not linearized), 
respectively. For all the calculations in this paper, we chose fre- 
quencies in the continuum with > 10'^ s~' to be linearized, 
while the rest were treated within the ALl scheme. Hence all 
the frequencies where Compton scattering could be important 
were linearized. In the future, we intend to implement an im- 
proved ALl treatment of Comptonization that will allow us to 
treat all frequencies as ALl frequencies, which would lead to a 
significant reduction of computer time. Work on this problem 
is in progress and will be described in a future paper. 

In the linearized transfer equation, there is no problem treat- 
ing the frequency coupling induced by Comptonization (i.e., 
the presence of the A, and C, terms in equation |l2|) or the 
non-linearity due to stimulated Compton scattering. However, 
before entering a linearization step, one must perform a "for- 
mal solution" of the transfer equation, i.e., a solution with the 
known, current thermal source function. This is very easy in 
the case of Thomson scattering because one can simply solve 
the transfer equation for one frequency at a time. The angu- 
lar dependence of the radiation field is also easy to deal with 
using the Feautrier method (e.g., Mihalas 1978). The formal 
solution determines the variable Eddington factor, fy = Ky/Jy, 
where Ky is the second moment of the specific intensity: Ky = 

{\|2)f_Jy{^i)d^i. 

In the case of Compton scattering, even the formal solution 
presents numerical complications because we cannot use the 
frequency-by-frequency strategy. In the present version of our 
computer program TLUSDISK (see also Paper 11), we use the 
following procedure for the formal solution of the transfer equa- 
tion with Comptonization. We first solve the transfer equation 
in the second-order form (e.g., Mihalas 1978), viz. 



dr} 



(19) 



where we assume a current value of the Eddington factor /, . We 
then discretize in depth, and organize the mean intensities 7, ^ 
(/ is the frequency index and d— 1 , . . . ,D is the depth index) in 
a set of column vectors 



Ji = {Ji.i,Ji.2, 



(20) 



so that the mean intensity vector contains intensities at all depth 
points for a given frequency point /. The resulting discretized 
transfer equation reads 



- A/J/_i + B,J/ — C/J/+1 — L; 



(21) 
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Here the matrices B are tridiagonal (because of the differ- 
ence representations of the second derivative with respect to 
depth), while the matrices A and C are diagonal (because the 
terms containing the frequency derivatives are local in physical 
space). Detailed expressions for the elements of matrices A, B 
and C are given in Appendix A. We avoid the nonlinearity due 
to the stimulated scattering term by writing the Compton source 
function as 

^comp, +BiJi+CiJi+i +Ji{UiJ°L'i +y,-/°'^ + , 

(22) 

i.e., by assuming that the stimulated contribution is evaluated 
using the current ("old") values of the mean intensity. If 
needed, the formal solution may be iterated to obtain an im- 
proved stimulated contribution, or one may devise a lineariza- 
tion procedure for the formal solution itself (which is distinct 
from our global linearization where we solve simultaneously all 
the structural equations). We have tested both approaches, and 
found that such refinements are never needed in actual model 
calculations. 

Solution of equation ( plj ) is done by standard Gauss-Jordan 
elimination, consisting of a forward elimination followed by a 
back-substitution, just as in the original Feautrier method for 
the Thomson scattering case. 

Finally, the new Eddington factors are obtained by perform- 
ing a set of elementary formal solutions for one frequency and 
angle at a time, with the total source function evaluated using 
the previously determined values of the mean intensity at all 
frequencies. 

2.2. Treatment of Metals 

We consider the following chemical species: H, He, C, N, O, 
Ne, Mg, Si, S, Ar, Ca, Fe, and Ni. We allow for all possible ion- 
ization stages, from neutrals to fully stripped atoms (e.g., Fe I 
to Fe XXVII). For the computations presented in this paper, we 
treat all ions except H I and He II as one-level atoms. This 
is done in order to explore the basic effects of metal opacities, 
together with Comptonization and departures from LTE, in an 
efficient way. We plan to use more extended model atoms in fu- 
ture calculations; very likely by adopting model atoms similar 
to those in the X-ray photoionization code XSTAR (Kallman 
2000). For hydrogen, we use a 9-level model atom (where the 
ninth level is a merged "superlevel" of all higher states) as in 
Papers I-III. He II is represented as a 14-level atom. Since we 
consider only hot disks here. He I is represented by a 1 -level 
atom, in contrast to Papers I-III where it was represented by a 
14-level atom. We tested that this approximation does not lead 
to any errors since He I has negligible abundance in the current 
models. 

For hydrogenic ions, we use the standard analytic expression 
for the photoionization cross-section. For all nonhydrogenic 
metal ion photoionization cross-sections, we use data from XS- 
TAR, version 1, kindly supplied by Tim Kallman (subroutine 
bkhsgo, based on data from Barfield, Koontz, & Huebner 
1972). XSTAR was recently upgraded to version 2, and now 
uses more recent sources of photoionization cross-sections (the 
Opacity Project and others); we also plan to adopt these cross- 
sections in future calculations. 

We take into account all inner-shell photoionization pro- 
cesses, but make the simplifying assumption that if an Auger 
electron is energetically possible, then it is in fact produced and 
the photoionization results in a jump by two stages of ioniza- 
tion to a ground state configuration. We therefore neglect both 



fluorescence and multiple Auger electron ejection arising from 
inner shell photoionization. Note that this is not a good ap- 
proximation for low ionization stages of metals, where multiple 
Auger electron ejection is likely after an inner shell photoion- 
ization (e.g. Kaastra & Mewe 1993). This might be important 
in regions where non-LTE effects are important in determining 
the abundance of these ion species. In future calculations we 
will lift these simplifying assumptions and consider more elab- 
orate inner shell photoionization branching possibilities. 

We include the effects of collisional ionization of the outer- 
most valence shell electron of each metal ion species by again 
using data from XSTAR, version 1 (subroutine cion). Once 
again, we plan on improving this treatment with more modern 
data in the future. 

Dielectronic recombination is handled using the subroutine 
direc from XSTAR, version 1. This is based on data from 
Aldrovandi & Pequignot (1973), Nussbaumer & Storey (1983), 
and Arnaud & Raymond (1992). Numerically, we treat dielec- 
tronic recombination by introducing an artificial modification 
of the photo-ionization cross-section, which has the advantage 
that the detailed balance relations are automatically satisfied. 
We describe this procedure in more detail in Appendix B. 

3. UNDERSTANDING TEMPERATURE STRUCTURE IN THE 
PRESENCE OF COMPTONIZATION 

Before we turn to actual numerical results, it is worthwhile 
developing a simple analytic model that will allow us to under- 
stand the vertical temperature profile in the presence of Comp- 
ton scattering. A similar study for the case of "classical" disks 
without Comptonization was done by Hubeny (1990). The ap- 
proach developed here is its generalization. We note that Madej 
(1989) considered a similar LTE-gray model in the context of 
neutron star atmospheres (i.e., with no viscous heating), al- 
though he only provided a numerical solution. 

3.1. Energy Balance 

The energy balance equation may be written in the form (see 
Paper II) 



dFnd 
dm 



eff " 



w{m) 
mow 



(23) 



Here w{m) is the dissipation rate at column density m and w 
is the vertical mean of w. In evaluating w we employ the rel- 
ativistic corrections of Riffert & Herold (1995) as in Paper II; 
replacing them with the forms advocated by Abramowicz et al. 
(1997) produces no significant change in the emergent spec- 
trum. 

The important point is that we can express the gradient of the 
radiation flux through the transfer equation 



rad 



dm 



4-Tr 



Xu{Sf~J,)dv. 



(24) 



Equating now the right-hand sides of equations (^) and ( p4| ) 
gives another form of the energy balance equation, which ex- 
presses the fact that the total energy generated locally is bal- 
anced by the net local radiative loss. 

The integrand of the radiative loss term is given by 



xJ^ + {x-3Q)Jl + QJ'^ + 



2hv- 



■Ju2x{Jl-Jy) 



.(25) 
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In the following, we neglect the stimulated scattering term, and 
assume that the thermal source function is given by the Planck 
function. We stress that this is done in this section only for de- 
veloping a simple analytic model; in actual calculations these 
simplifications are not made. The integral in Eq. ( |24l ) can be 
evaluated using integration by parts for the derivative terms. Af- 
ter some algebra we obtain 



kbB - KjJ + K„ (49 -x)J=^ 7;tf ^ 

47r mow 



(26) 



where J = Ji,dv, and analogously for B. The mean opaci- 
ties are defined by 



Jo 



JQ 

(which is called the Planck-mean opacity), and 

Kcr = n^CFll p. 



Finally, x is defined by 



h _ h vJjjdv 



(27) 



(28) 



(29) 



(30) 



which has the meaning of a mean frequency weighted by the 
local mean intensity. Since the mean intensity is typically a 
sharply-peaked function of frequency, D is close to the fre- 
quency where J^, attains its maximum value. We note that the 
traditional Compton temperature, 7c, is defined (e.g., Krolik 
1999) as the temperature for which the term 40 — x vanishes, 
i.e. 

hv 



Ak 



1.2 X 10" 



(31) 



Equation ( [26[ ) has a simple physical interpretation. The right- 
hand side expresses the net energy generated per unit mass at a 
given depth in the disk. The left-hand side expresses the net en- 
ergy loss, i.e. the energy that has been transferred from material 
particles to photons. The term kbB — hjJ expresses the trans- 
fer of energy from particles to photons by thermal processes 
(i.e. bound-bound, bound-free, and free-free atomic processes), 
while the term (40 — x)J expresses the net loss of energy of 
particles, in this case electrons, by Compton scattering. The 
latter term is in fact a generalization of the usual formula (e.g., 
Rybicki & Lightman 1979) for the net photon energy gain (i.e., 
net electron energy loss) in a single scattering, which in our 
notation reads 

Mhv) 

— ^ — -=A<d-x. (32) 
hv 

We also note that the ratio Ks/na can be expressed as 

POO />oo 

HB/i^cy = I {K^^/{neaj)B^diy/BK / e^Bi,dv/B = e, 
Jo Jo 

(33) 

i.e. the ratio of the Planck-mean to the scattering opacity has the 
meaning of a mean photon destruction parameter. Notice that 
the last approximate equality only holds in the case of dominant 
electron scattering, e^<^l. 



To proceed further, we need an expression for the local inte- 
grated mean intensity, J. This can be obtained from the second 
moment of the radiative transfer equation 



dK^ 
dm 



(34) 



where H^, — (1/2) J^^I^{p.)dp is the first moment of the spe- 
cific intensity. We now invoke the Eddington approximation, 
i.e., Ki, = and integrate equation ( p4| ) over all frequencies. 
We obtain, again in the case of dominant electron scattering, 

dJ , „ 3._ ^3^^ 



dm 



4-TT 



The local radiation flux follows from integrating equation (p3|), 
viz. 

w{m')dm' 



1-- 



mow 



(36) 



Substituting equation (^6|) into ( |35| ) and solving it for J we ob- 
tain 



(37) 



where we used the Eddington approximation form of the 
boundary condition, 7(0) — H{0)\/3, and expressed H{0) 
through the effective temperature. The optical depth is defined 
by 



Kadm' W Ka^m 



(38) 



since we assume that hydrogen and helium are completely ion- 
ized so that Kcr is nearly constant. This optical depth is thus 
the Thomson optical depth. The other optical depth in equation 
( ^7| ) may be called the "viscosity-weighted" optical depth, and 
is defined by 



dm' K„ I w{m")dm" /mow . (39) 
^0 

In the case of a depth-independent kinematic viscosity. 



T 

2rto 



(40) 



where Ttot = K^-mo is the Thomson optical depth at the mid- 
plane. 

Finally, combining equations ( ^ ) and (^7|), we arrive at the 
basic equation of this section, 

f T 3 f I \r M w(m) 

^'[l^J -4[v5^^-V^^'-^'^'''-''^^''^o^ 

(41) 

where we used the well-known expression B = {a/Ti)T'^. 

Using equation ( 03|), and assuming, as usual, that Ky = Kg, 
we rewrite equation ([4l|) in a more instructive form. 



'eff 



1 



4 VV3 



T-r„, (£-46 



w{m) 1 

W AT,r 



(42) 

We see that the effects of Compton scattering are negligible 
if |40— x| ^ e, in which case we recover equation 3.11 of 
Hubeny (1990). 

Before we discuss the meaning and implications of equation 
( p^ in more detail, we first turn to estimating the mean fre- 
quency. 
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3.2. The Mean Frequency 

A determination of the mean frequency is easy if the mean in- 
tensity is proportional to the Planck function, J^, ^ aB,y, where 
a is a frequency-independent constant. In this case 

- _ /o°° ^V[sxp{h,y/kT) -\]dy ^ kT r(5) C(5) kT 

j^vy{ti^p{hv/kT)-\]dv /i r(4) c(4) ~ ■ h' 

(43) 

Here, F is the Euler gamma function and C is the Riemann 
zeta function (e.g., Abramowitz & Stegun 1970). The Compton 
cooling term is thus given by 

4e-jc = e(4-3.83) = 0.17e = 2.867x 10""r. (44) 

Deep in the atmosphere, 7^ w B^{T), where T is the local tem- 
perature. Therefore, the Compton cooling term is positive, and 
is determined by the local temperature. 

At the surface, the mean frequency can be estimated as fol- 
lows. Let us assume for simplicity that the total Thomson op- 
tical depth is large. Then the surface mean intensity, assuming 
the Eddington approximation, is given by (see, e.g., Mihalas 
1978 or Rybicki & Lightman 1979) 

^(r = 0)« V^B.(r(T=l/V3^)). (45) 

In other words, the mean intensity is still proportional to the 
Planck function; however the constant of proportionality, -/eU, 
depends on frequency, and the temperature at which the Planck 
function is evaluated is not the local temperature, but the tem- 
perature at depth r = r^°™ = 1 /v^3e^, which may actually be 
quite large. We call the depth r^"™ the effective depth of for- 
mation. We note that there is an intimate relation between the 
effective depth of formation and the traditional effective optical 
depth (e.g., Rybicki & Lightman 1979; their equation 1.98). 
The effective optical depth is given by (in our notation) 

< - \fr^u ~ V^- (46) 

Our effective depth of formation is thus equal, up to a factor 
of the order of unity, to the optical depth of the point where 
the effective optical depth is equal to unity, i.e., r* w 1 at 
Ti/ ~ 1 1 \ffZ- In view of all the approximations made in deriving 
Eq. (Esf) (the Eddington approximation, a depth-independent 
ejy, etc.), the factor \fl> is inconsequential. 

To proceed further, let us assume that in the deep layers we 
have \Q~x — l.Zdl x 10^"r < e, and thus the Compton term 
can be neglected. This is actually a reasonable approximation 
for estimating the surface value of v for our disks, since typ- 
ical local temperatures are of the order of a few times 10'' K, 
while e « lO"'* in the deep layers (note that in the deep layers 
e is larger than near the surface because of the higher density). 
Assuming also Ttot 3> 1, we are left with a rough estimate of 
the temperature in the deep layers, T^ij) « (3/4)7'gffT, which 
follows from equation (|4l|), neglecting here for simplicity a cor- 
rection term of the order of unity coming from r,i . 

We denote the temperature at which the Planck function in 
equation ( p3[ ) is evaluated as T*; it is thus given by 

7^* = (3/4)i/47;jj3-i/8g-i/8^ (47) 
The dependence of T* on frequency arises only through the 

— 1/8 

term ty ' , and is therefore rather weak. We may then replace 



e^j'^^ by (e)^'/^. Similarly, if we assume that y/eZ varies with 
frequency much more slowly than fi^ (at least around the max- 
imum of the Planck function that contributes most to the inte- 
grals), then we may again use equation (^3]) to obtain 

P = 3.83 {kT*/h) = 3.83 (k/h) (3/4)'/'*3-'/** {ly'^^ns 

w6.5 X 10"'(e)"'/'*reff. (48) 

The Compton cooling term may thus be expressed as 

k 

46 -x= ^ (4r- 3.83 r*) =6.75 x 10""'(r -0.96r*) , 

nigC-^ 

(49) 

where we left the "formation temperature", T*, unspecified. 
This equation is instructive because we easily see that if T > 
0.96 T*, Compton scattering leads to cooling of the surface lay- 
ers, while if r < 0.96 T*, Compton scattering heats the layers. 
If Compton scattering is the dominant source of opacity and 
emissivity, then the surface layers have the temperature T k,T* . 

3.3. Surface temperature 

An estimate of the surface temperature is of considerable in- 
terest since it will show whether our models are in principle 
able to produce high-temperature external layers, traditionally 
called disk coronae. 

Using the results of the previous section, we can rewrite 
equation ( p^ at the surface (i.e. r = t^. = 0) as 

f TqV V3 r 4;t , ^ w(0) 1 

At^] =^ ' j(r„-o.96r) +^^. 

(50) 

Using equation (p^, and introducing a normalized surface tem- 
perature, t = T()/ Teff , we obtain the following fourth-order alge- 
braic equation for t, 

at^ + bt = c, (51) 

where 

a = e, ^ = 2.92x lO^'^Jeff, (52) 

and 

c = 0.43e + 2.92x 10-'°7;ff (e)"'/** + — . (53) 

w 4rtot 

A general solution of equation ( |5l] ) is possible to find analyt- 
ically, but is quite cumbersome. However, we may consider 
several limiting cases: 

l)fl>Z7,i.e. e>3x lO^'^Teff 
This is the limit of negligible Compton scattering, which ap- 
plies if either the mean photon destruction parameter e is very 
large, or the effective temperature is very low. In this case the 
solution of equation ( ^1| ) is simply t — (c/a)'/"^, i.e. 

To f w{Q) 1 

0.43 + ^ — ^ . (54) 

There are two limiting cases in this situation. If the second term 
in the brackets in negligible, then one is left with the classi- 
cal stellar-atmosphere surface temperature: 7b ~ 0.81 T^sf. The 
second term is important, if (i) w{Q)/w is large, i.e., if the vis- 
cosity increases towards the surface; or (ii) if the term ertot 
is small. Note that the latter term is in fact the total Planck- 
mean optical depth, i.e. the total optical depth for pure thermal 
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absorption. Therefore, the surface temperature may be large 
(i.e., larger than T^ff) even in the absence of Compton scatter- 
ing and an outward-increasing dissipation. This is readily ex- 
plained. Thomson scattering does not influence the energy bal- 
ance. Since we now have an optically thin disk with respect to 
the total Planck-mean opacity, the total energy radiated away is 
proportional to St, where S is the source function, and t the 
total optical depth. We take S = B (x T'^, while the total energy 
dissipated is proportional to r^"^. These energies are equal, so 
that we obtain {T /T^ff)* oc l/r. In other words, an effectively 
thin slab must be hotter than a thicker slab to be able to radiate 
the same dissipated energy. 

2)fl<Z7,i.e. e<3x lO-'^Teff. 
This is the limit of dominant Compton scattering. The solution 
of equation dsTj ) is f = c/b, i.e. 

provided that at'* is also <C bt, i.e., t <^ {b/aY^^. When this 
criterion is satisfied, there are two limiting cases. If the second 
term on the right-hand side is small, then Tq w (e)^'/** T^ff. For 
typical values of e of the order of lO"*" to 10"'', we see that 
Compton scattering leads to a heating of the surface layers up 
to a factor of 3 - 6 above the effective temperature! Obviously, 
the surface temperature is larger for smaller values of e because 
Compton scattering becomes more efficient with respect to the 
thermal processes. On the other hand, if the second term is im- 
portant, either because w(0)/vi> is large or the total optical depth 
is small, then the solution is given by the second term of equa- 
tion (|55|). When of"* is not small compared to bt, equation ^54| ) 
is again roughly correct. 

We summarize the behavior of the surface temperature in 
Figures 1 and 2. 




Fig. 1- Approximate surface temperature, expressed as To/Tcff, dis- 
played as a function of T^ft, for different values of e, and for two values 
of Tiot. The individual curves are labeled by the value of log e. Solid 
lines are for a total optical depth rto, = 10*, and dotted lines are for 
r.o, = 10\ 

For simplicity, we assume constant viscosity, i.e. w{0)/w = 
1 ; we will study the effects of an increasing viscous dissipation 
towards the surface in a subsequent paper Figure 1 displays the 
surface temperature, expressed as Tb/Teff as a function of Teff, 
for different values of e, and for two values of Ttof Figure 2dis- 
plays the surface temperature as a function of the total optical 
depth, for several values of e and Jeff- 
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Fig. 2- Approximate surface temperature, expressed as To/Tcff, dis- 
played as a function of rtoi, for different values of e, and for three val- 
ues of Tcff. The individual curves are labeled by the value of e and by 
the value of logTcff. Solid lines have Tcff — 10^ K, dotted lines have 
Jeff = 10^ K, and dashed lines have Teff = 10'^ K. 

We now turn to the accurate numerical models, and show that 
the above estimates agree reasonably well with exact numerical 
results. 



4. RESULTS 

4. 1 . Comparison to the models from Paper III 

We shall first present a comparison of models computed with 
and without Comptonization. We take the hottest disk from our 
grid of Paper 111, where the effects of Comptonization are ex- 
pected to be the largest. This is a disk around a black hole with 
M — {\ /8) X 10'' M0 and mass accretion rate M = 1/4 Mo/yr 
(i.e. with luminosity L « 0.3LEdd)- To isolate the effects of 
Comptonization, we assume, in agreement with Paper III, that 
the disk consists of hydrogen and helium only. Also following 
Paper 111, we take the local viscosity to be described by a two- 
step power law, with a constant kinematic viscosity in the inner 
99% of the disk mass, while the viscosity is assumed to de- 
crease as m^^/-' in the outermost 1% of the mass (see Paper 111 
for details). The reason we chose an artificial decrease of vis- 
cosity with height was to prevent a "thermal catastrophe" of the 
disk in the low optical depth regions (e.g., Shaviv & Wehrse 
1986). In the case of H-He models without Comptonization, 
there are no efficient mechanisms that are able to radiate the 
dissipated energy away when the atmosphere becomes optically 
thin in the most opaque transition; the atmosphere thus under- 
goes a thermal runaway. This is obviously prevented when one 
considers Compton instead of Thomson electron scattering, as 
we have shown analytically in the previous section, and will 
demonstrate numerically in this and subsequent sections. We 
show the models with vertically decreasing dissipation only for 
the purpose of direct comparison of the present results to those 
of Paper 111. 
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Fig. 3- Temperature as a function of depth for selected annuli. From 
top to bottom, the curves correspond to radii R/Rg — 2, 3, 4, 5, 6, 7, 8, 
9, 10, 12, 14, 16, 18, 20, 25, 30, 40, 50, 60, 70, 80, and 90. Thin lines 
represent the models from Paper III computed without Comptoniza- 
tion, while thick lines correspond to models with Comptonization. 

Figure 3 displays the local electron temperature as a function 
of position for selected annuli. The position is expressed as 
column mass, m, above a given depth. The models with Comp- 
tonization exhibit a temperature rise at the surface, even for the 
case of viscosity decreasing outward. This is easily explained 
by the considerations presented in the previous section. The 
hottest annuli have T^ff « 2 x 10^ K. The typical value of e is 
about 10"^ so that e < 3 x lO^^'Teff « 6 x 10"^ and thus 
we are in the regime of important Compton scattering, where 
Eq. (55) applies. Since viscosity is forced to decrease outward, 
the second term on the right-hand side of Eq. ( |5^ is negligi- 
ble, and we therefore find T /T^ff w 2, which is indeed seen in 
Figure 3. The cooler annuli have progressively larger values 
of e and lower T^ff, and therefore are less and less influenced 
by Comptonization. Again, this is nicely verified by numerical 
calculations. 
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Fig. 4- Temperature as a function of depth for selected individual an- 
nuli. The curves correspond to radii R/Rg (from top to bottom) 2, 3, 
4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 25,' 30, 40, 50, 60, 70, 80, and 
90. Thin lines represent the models with constant dissipation, while 
thick lines represent the models with decreasing dissipation towards 
the surface. 

To see the effect of changing the dissipation law on the 
temperature structure, we have computed the same models as 



above, but keeping the kinematic viscosity constant with depth, 
i.e. without an artificial decrease of viscosity at the surface. 
The temperature structure for the individual annuli is displayed 
in Figure 4. Note that models with constant dissipation were 
impossible to construct with the original assumption of Thom- 
son scattering, because the temperature exhibited a runaway 
instability behavior. The behavior of the surface temperature 
is again easily explained using Eq. (55). In this case, we 
have w(0)/vT> — 1. Moreover, for the hottest annuli we have 
Teff w 2 X 10^ K, and « 2 x 10^, so the second term in equa- 
tion (55) is about 2. Since the first term is also roughly equal 
to 2, we see that in this particular case the Compton heating 
and the viscous heating have a comparable effect on the surface 
temperature. The resulting T /T^ff should be around 4, which is 
again nicely verified by numerical calculations. 
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Fig. 5- Local emergent flux for the individual annuli. From top to 
bottom, the curves correspond to radii R/Rg — 1.5, 2, 2.5, 3, 3.5, 4, 5, 
6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100, 
120, 140, 160, 180, and 200. Thin lines represent the models from 
Paper III computed without Comptonization, and thick lines represent 
the models with Comptonization. 

From a practical point of view, the most important feature 
is the behavior of the emergent flux from the individual an- 
nuli when Comptonization is taken into account. Figure 5 dis- 
plays a comparison between the model without (thin lines) and 
with (heavy lines) Comptonization for the individual annuli dis- 
played in previous figures. 

The models with constant and with decreasing kinematic vis- 
cosity at the surface are indistinguishable in the plot. This is 
because only the surface layers are influenced by changes in 
the dissipation law, while the emergent continuum radiation is 
formed in much deeper layers. The emergent flux would be 
influenced if the dissipation increases outward, in which case 
a hot corona (with temperatures of the order 10^ to 10^ K) is 
formed. Such a corona would produce a significantly increased 
X-ray flux. We will present a more extended study of this phe- 
nomenon in the next paper of this series. 

Comptonization in hot annuli leads to the well-known in- 
crease of the flux in the EUV and the soft X-ray region (v < 
10'^ s^'). For larger frequencies, the flux drops fast, so that the 
flux for the Comptonized model is lower than the flux for the 
original, Thomson-scattering, model. An analogous behavior 
was also found for the emergent spectra of (hotter and denser) 
neutron star atmospheres (Madej 1991). This behavior is read- 
ily explained by examining the net photon energy gain (or loss) 
as a function of frequency. As discussed in Section 3 (see also 
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Rybicki & Lightman 1979), the net energy gain of a photon 
with frequency v is given by 40 — x. Since we are dealing with 
a depth-dependent situation, one has to evaluate the appropri- 



ate 8 at the effective depth of formation, Tjj w t, 



form , 



1/ 



Note that when r*°™ > rtot, where rtot is the optical depth 
at the midplane (generally, it is a frequency-dependent quan- 
tity, but the dependence is very weak in the case of dominant 
electron scattering), the disk is effectively thin. The effective 
depth of formation should thus be more generally given by 
^foim ^ min(l/y^, Ttot). Nevertheless, we will assume an ef- 
fectively thick disk in the following discussion. The tempera- 
ture at the formation depth is thus given by [see also Eq. (47|)], 



' Teff El/ 
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Fig. 6- Computed value of 40(7^*) compared to x as a function of 
frequency. The Compton parameter jcff is also displayed. All the quan- 
tities are dimensionless numbers. 

In Figure 6 we display the computed value of 48 (rj') com- 
pared to X as a function of frequency. We see that for v ^ 
6 X 10'*, the term 48 — x is positive, and therefore Comp- 
tonization leads to a net increase in photon frequency, while for 

^ 6 X 10'* the photon frequency is decreasedhy Comptoniza- 
tion. These values tell us how the photon frequency is shifted, 
but do not tell us about the number of available scatterings and 
thus about the resulting effects of the photon frequency shifts 
on the emergent spectrum. In order to verify that Comptoniza- 
tion in fact has a significant effect on the spectrum, we also 
show the effective Compton parameter yeff , which is defined by 
(Rybicki & Lightman 1979) 



3;eff = 48max[r;°™,(rf™)^] 



(56) 



Here the factor max[rfo„n, (Tfoim)^] expresses the total number 
of scatterings of a photon with frequency i/. Since in our case 



^form , 



1 /y/eZ and thus r*°™ > 1, we have 



3'eff«4x lO-'^Teffe 



-9/8 



(57) 



We can roughly estimate and therefore y^-g as follows. Let 
us assume for simplicity a relatively hot plasma where all hy- 
drogen and helium are almost completely ionized, and that the 
only sources of thermal opacity are the Lyman continua and 
the free-free transitions of hydrogen and ionized helium. The 
bound-free and free-free opacities are thus given by (see, e.g., 
Mihalas 1978) 



1.17 X \0^*T-^^^n,i^ 



-3 



0.158/r6 



H{iy,iyn) 



74 „ „o.m/T6 

L nHellie 



(58) 



and 



Kff(i/) w3.69x lO^r-'/^Hei'-^K+Z^nHeiii], (59) 

where Z = 2 is the effective charge of the ionized helium; Wp 
and «Heiii are the number density of protons and of helium 
nuclei, respectively; and = T/IQ^. Here we assume LTE 
for the hydrogen and ionized helium ground-state populations. 
We also assume that the bound-free and free-free Gaunt fac- 
tors are equal to unity. H{v,v\^) is the unit step function, de- 
fined in such a way that H = for ly < vq, and H = 1 for 
1/ > i^o; i^u = 3.29 X 10'^ s~' is the frequency of the hydro- 
gen Lyman edge. Since we assume almost completely ionized 
plasma, then from the charge conservation equation we obtain 
Hp = ne/(l +2y) and nHeiii = rif^Y /{I + 2Y). Substituting these 
values, and assuming the solar abundance of helium, F = 0.1, 
we obtain for 

w6.5xl0 nuT. i/,, 7^.(7) , (60) 

where nu = "e /10'\ z^i6 = i//10'*, and 7 is a function of order 
of unity that depends only weakly on T and v (at least in the 

temperature range T ^2x 10^ K relevant here). 
0.226 r 



7.(r) = l- 



T6 



(61) 

Equation (^0|) is valid for any temperature and electron den- 
sity. We may also specify the corresponding e at the effec- 
tive depth of formation, e* = €^{T*), where T* is given by 
T* ~ 7;ff(e*)~'/^ Substituting for e from equation (^, we 
obtain 

(62) 



r* ~3fiT'*/'5„-2/15 2/5 



[notice that the scaling 8(7'^') cx v'^^^ is nicely seen in Figure 6], 
and thus 



f* ~3 2X 10-57.-8/15 16/15 -16/5 



(63) 



Here we use the same notations as before, namely Teff e 



Wio*, and r;^ = r;/io' 



have 



Finally, from equation (q7|) we 

(64) 



-44x24/15 -18/15 18/5 



Again, the approximate scaling jeff oc ly^^/^ is verified in Fig- 
ure 6. 

We see from Figure 6 that jeff ~ 1 for ly ^ 3 x 10'*, and 
thus Comptonization should significantly influence the emer- 
gent spectral energy distribution for these frequencies. From 
Figure 6 we deduce that the emergent flux for frequencies 

3 X 10'* ^ ly ^ 6 X 10'* should be enhanced over the model 
with Thomson scattering only, while the opposite should ap- 
ply for higher frequencies, which is indeed demonstrated by 
the exact results shown in Figure 5. In fact, the enhancement 
of flux caused by Compton scattering is seen already for lower 
frequencies, which is a consequence of the higher temperature 
for the Comptonized model (cf. Figure 3). 

An inspection of Figure 5 shows that the emergent flux is 
influenced only for the hottest annuli, those with effective tem- 
perature above 150,000 K; the cooler annuh are influenced only 
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marginally. This is easily seen from equation (p4|). The ef- 
fective temperature varies with radial coordinate roughly as 



Teff oc R^^/'^, so that T^^^ cx R^^/^ . Estimating the appropri- 
ate electron density is much more complex. We consider two 
limits. In completely radiation-pressure dominated disks, and 
neglecting relativistic corrections, we may consider the mass 
density at the effective depth of formation to be equal to the 
vertically averaged mass density, po- Thus oc po oc R^^^. In 
this limit, jeft oc R^^, i.e., it decreases sharply with increasing 
radial distance. The opposite limit is to assume the electron 
density roughly constant or only slightly increasing with R (see, 
e.g. Figure 4 of Paper 111), which follows from an interplay be- 
tween the higher midplane density and a sharper decrease of 
vertical density profile for increasing radial coordinate. Even 
in this extreme case, yeff decreases faster than linearly with in- 
creasing R. These considerations explain why the Comptoniza- 
tion influences mostly the inner annuli, while the outer annuli 
are influenced progressively less and less. 

To see the effects of Comptonization on the hydrogen and 
He 11 Lyman edges more clearly, we present in Figure 7 a blow- 
up of Figure 5. It is clear Comptonization indeed considerably 
smears the He 11 Lyman edge. The hydrogen Lyman edge is 
smeared as well, although the effect is much smaller because 
Comptonization is much less efficient at the relatively low fre- 
quencies corresponding to the hydrogen Lyman edge. 
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Fig. 7- A blow-up of Figure Sshowing the region of the hydrogen and 
He II Lyman edges. From top to bottom, the curves correspond to local 
emergent flux for the individual annuli at radii R/Rg — 2, 3, 4, 5, 6, 7, 
8, 9, 10, 12, 14, 16, 18, 20, 25, and 30 Thin lines represent the mod- 
els from Paper III computed without Comptonization, and thick lines 
represent the models with Comptonization. 

Finally, Figure 8 presents disk-integrated spectra for the 
above disk model. As in Paper 111, the spectrum is found by in- 
tegrating the total emergent intensity over the disk surface using 
our relativistic transfer function code (Agol 1997). The trans- 
fer function computes the trajectories of photons from infinity 
to the disk plane, finding the radius, redshift, and intensity at 
each angle at infinity for a given observation angle (Cunning- 
ham 1975). We neglect the effects of radiation that returns to 
the accretion disk. 
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Fig. 8- Integrated spectral energy distribution i/L^ (in erg s ') for a 
disk withM = (1/8) x lO' Mq, M = 1 /4M0 yr~', (i.e. L/LEdd = 0.3), 
and cos / = 0.5, for the model with (bold line) and without (thin line) 
Comptonization. 

Comptonization hardly affects the optical and ultraviolet por- 
tions of the integrated spectrum, but it has a more dramatic ef- 
fect at higher energies. The flux distribution around the He 11 
Lyman limit is generally smoother; the flux in the region be- 
tween 2 X 10"' and 5 x 10"' Hz (i.e., roughly between 0.1 and 
0.2 keV) is about 30% larger if Comptonization is taken into 
account. About 10 percent of the total luminosity is shifted 
from the high energy tail (above 6 x 10'^ Hz) to the peak of the 
spectrum (between 1.5 and 4 x 10"' Hz.) 

We stress that, because they had larger central masses and 
therefore lower central temperatures, the effects of Comp- 
tonization for all other disks presented in Paper III should be 
smaller Therefore, we conclude that neglecting Comptoniza- 
tion made little difference to the grid of models presented in 
that paper. This conclusion is in general agreement with the re- 
sults of Laor & Netzer (1989), who found Comptonization to 
be unimportant for disks with T^s ~ 10^ K. 

4.2. Comparison with Previous studies 

Studies of the influence of Comptonization on disk atmo- 
spheres were pioneered by Ross, Fabian, & Mineshige (1992 
- hereafter referred to as REM). We now present a comparison 
with the models of REM, who computed models of AGN accre- 
tion disks taking into account the effects of Compton scattering 
and some NLTE effects. Specifically, they assumed a constant 
(depth-independent) density in the z-direction, and solved the 
energy balance equation, the radiative transfer equation (but in 
the diffusion approximation), and the rate equations for the first 
three levels of H 1 and He II. In the latter equations, the radia- 
tive rates in the bound-bound transitions were treated by means 
of escape probability theory. 

In order to compare our treatment of Comptonization to 
theirs, we have computed a set of test models as similar to REM 
as possible. We have kept density fixed, and consider only three 
levels for H I and He II. Otherwise we have used our procedure 
to solve all the structural equations except hydrostatic equilib- 
rium. Since our program does not use escape probabilities, we 
have determined the bound-bound radiative rates exactly, by 
solving the transfer equation with 15 frequency points per line. 

We compare a single annulus at R/Rg = 9.8, taken from 
a disk with M = 10^ Mq (Schwarzschild black hole) and 



11 



i/iedd = 0.3 (i.e. M ^O.l Mo/yr) to the corresponding RFM 
model, kindly supplied to us by Randy Ross - see Figure 9. 
Results for other annuli are analogous. The agreement of the 
computed flux in the continuum is very good, although there 
are some differences, particularly at the shortest and the longest 
wavelengths. The difference in the region of the flux peak is 
about 10%, which is likely explained by the RFM treatment of 
the inner boundary condition for the radiative transfer equation. 
For the frequencies at which the total effective optical thickness 
at the disk midplane is larger than unity, RFM assume that the 
radiation at the effective optical depth equal to unity is in LTE 
at the temperature given by an approximate expression [their 
Eqs. (15) and (16)]. In the present case, the effective optical 
thickness is slightly above unity for u < 10^^, and reaches val- 
ues of about 10 at the He II Lyman discontinuity. Therefore the 
RFM treatment of the inner boundary condition is somewhat 
inaccurate for these frequencies. 

The predicted profiles for the hydrogen and helium lines dif- 
fer significantly, and some are even flipped from emission to ab- 
sorption. This is very likely explained by differences in treating 
the lines, in particular the escape probabiUty treatment adopted 
by RFM. Another reason may be a rather coarse wavelength 
resolution, and consequently low number of frequency points, 
adopted by RFM. As mentioned earlier, we will perform a de- 
tailed study of the effects of spectral lines on the vertical struc- 
ture and emergent spectra for AGN disks in a future paper of 
this series. 
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Fig. 9- A comparison of a sample model for one annulus at R/Rg = 
9.8, for a disk with M = lO' M© and L/LeM = 0.3 (i.e. M = 
0. 1 Mq I yr) computed by TLUSDISK (solid line) with the correspond- 
ing Ross, Fabian, Mineshige (1992) model (dotted line). The abscissa 
is the total flux from the annulus, i.e. the flux per unit area times the 
total area of the annulus, which is chosen to extend from RfRg = 9.2 
to 10.4. 

In a subsequent study, Shimura & Takahara (1993) devel- 
oped a more accurate treatment of Comptonization. They con- 
sidered, as in the present paper, the Kompaneets form of the 
Compton scattering emissivity, and solved the same list of 
structural equations as in this paper (hydrostatic equilibrium, 
energy balance, and the radiative transfer equation). How- 
ever, they considered only free-free opacity, ignoring the many 
atomic opacity mechanisms, and they assumed a two-stream 
approximation for the radiative transfer (Rybicki & Lightman 
1979). Unlike the present paper, they assumed separate energy 

^We used the same prescription in Paper in. Note that there were typos in i 
chosen to be 5000 K. 



balance equations for the protons and the electrons. However, 
they also showed that there was very little difference between 
the two temperatures, supporting our approach of having the 
same kinetic temperature for aU particles. 



4.3. Effects of Comptonization in hot disks 

We have constructed several disk models. We consider three 
representative black-hole masses, M = 10*", 10^, and 10** Mq. 
For all three masses we take L/Lndd ~ 0.3, i.e. mass accretion 
rates equal to 0.002, 0.02, and 0.2 Mq yr~S respectively. We 
assume a maximum rotation Kerr black hole. For each disk, we 
consider two values of the viscosity parameter a: 0. 1 and 0.01 . 
Analogously to Paper III, we take the individual annuli for the 
radial coordinates R/Rg (with Rg = GM/c^ being the gravita- 
tional radius) equal to 1.5, 2, 2.5, 3, 3.5, 4, 5, 6, 7, 8, 9, 20, 12, 
15, 20, 25, 30, 40, 50 ,60, 70, 80, 90, 100, 150, 200, 250, 300, 
400, 500, 600, 700, 800, 900, 1000, 1200, and 1500. A ver- 
tical structure model is computed for all annuli for which the 
effective temperature is higher than 5000 K; the cooler annuli 
are assumed to radiate as black bodies.^ The outer edge of the 
disk is chosen to be the radius at which Teff reaches 1000 K. 

We first show the temperature structure for four representa- 
tive annuli of the disk with M = 10*" Mq, M = 0.002 Mq yr-\ 
and a = 0.1 and 0.01 - see Figure 10. 
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Fig. 10- Temperature as a function of depth for several representative 
annuli of flie disk for M = 10* M®, L/Lndd w 0.3; with a = 0.01 (soUd 
lines), and a = 0.1 (dashed lines). From top to bottom, the curves cor- 
respond to radii RjRg = 1.5, 3.5, 5, and 10. 

The figure illustrates several interesting features. First, it 
shows the dependence of the surface temperature on both ef- 
fective temperature and total optical thickness, and thus on the 
value of a. We recall that the total thickness of the disk is in- 
versely proportional to a (Paper III). The hottest annulus has an 
effective temperature w 8 x 10^ K. The total optical thickness 
at the disk midplane is approximately 240 for the a = 0.1 disk, 
and 2400 for the a = 0.01 disk. The mean photon destruction 
parameter for the a = 0.1 disk is about 10~^, while that for the 
a = 0.01 disk is about 10~^. From Figure 2 we deduce that 
To/Tgff should be around 6 and 2, respectively, for these two 
disks, which is indeed verified by numerical calculations. 

paper: the transition temperature was stated to be 4000 K, whereas in fact it was 
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the region around 2 x 10'^ Hz, so that the flux for lower fre- 
quencies is depressed compared to the model with a = 0.01. 
For higher frequencies, Comptonization is in contrast more ef- 
ficient for the model with a = 0.01, which leads to a larger shift 
of photon energies towards lower frequencies, thus resulting in 
a lower flux for the model with a = 0.01. 



Fig. 11- Local emergent flux for several representative annuli af the 
disk for M = 10*^ Mq, L/haid ~ 0.3; with a = 0.01 (solid lines), and 
Q = 0. 1 (dashed lines). The curves correspond to isidn R/Rg (from top 
to bottom) 1.5, 3.5, 5, and 10. 

Finally, the corresponding emergent flux for the four annuli 
displayed in Figure 10 are shown in Figure 11. Higher val- 
ues of a have somewhat suppressed flux near the peak and en- 
hanced flux at low and high frequencies. This effect is most 
pronounced at small radii, i.e. lower total optical thickness. To 
understand the behavior of the emergent flux, we display the 
effective Compton y-parameter and other interesting quantities 
for the annulus at R/Rg 1.5 of the disk with M = lO*" M©, 
i/LEdd ~ 0.3, and for the two values of a - see Figure 12. This 
plot is analogous to Figure 6. 
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Fig. 12- Computed values of 40(7^"), x, and the Compton parame- 
ter ycff for the annulus at R/Rg = 1.5 of the disk with M = 10* Mq, 
L/LEdd ~ 0.3, and for a = O.oi and 0.1. 

Increasing a decreases density in the disk, leading to greater 
Comptonization because jeff oc n^'^^'^ - see equation (|64|). In 
line with this expectation, jetf in the models with a = 0.1 is 
more than an order of magnitude greater than in the models 
with a = 0.01. This is indeed nicely seen in Figure 12. The 
frequency where 40(7'^*) w x is about 2 x 10'^ Hz. However, 
the disk with a = 0. 1 becomes effectively thin at i/ w 1 .3 x 10'^ 
Hz, while the annulus with a = 0.01 is effectively thick for all 
frequencies, so that jeff keeps increasing with frequency. In the 
model with a — 0.1, Comptonization is more efficient at re- 
moving photons from the frequency range < 2 x 10'^ Hz to 



4.4. Effects of metals 

We first computed a series of test models in order to deter- 
mine which metals must be treated in detail. These models had 
the full array of metals as discussed in Sect 2, i.e., H, He, C, N, 
O, Ne, Mg, Si, S, Ar, Ca, Fe, and Ni. These test models showed 
that Mg, Si, S, Ca, and Ni have only a marginal influence on 
models (both on the structure and emergent spectral energy dis- 
tribution), so we have computed most models with only H, He, 
C, N, O, Ne, and Fe. 

Figure 13 shows a comparison between models computed as- 
suming a simple H-He chemical composition and models with 
metals. We plot here the total integrated spectrum of the three 
representative disks. We also show the integrated spectrum as- 
suming a black-body flux for each annulus. As in our previous 
models, the self-consistent treatment of radiation leads to a re- 
distribution of flux from the optical and UV region to the EUV 
and and soft X-ray region as compared to the integrated black- 
body energy distribution. 

The effect of metals is relatively small; the only difference 
in the predicted spectra is seen in the high energy tails. Gen- 
erally, the spectra begin to differ at a frequency about three 
times higher than the frequency where the emergent flux attains 
its maximum. The reason for this may be seen more clearly 
in Figure 14, where we present a comparison of the effective 
Compton y-parameter for a representative H-He model and for 
a model with metals, for a representative hot annulus of the 
disk with M = lO*" Mq, L/Lndd » 0.3, and a = 0.01. We do not 
show a comparison of local temperature for the models with 
metals to those of the H-He composition because the differ- 
ences are very small. The metals thus do not influence the emer- 
gent radiation through their effects upon the vertical structure, 
but only through their effect on the thermal coupling parameter 
e and consequently on the efficiency of Comptonization. 

We see in Figure 14 that in contrast to the H-He model, the 
effective Compton y-parameter reaches a maximum in the re- 
gion 10i7~i/-2x 10'^ because of the effects of metals, in 
particular the O Vlll Lyman edge at 1/ « 2 X 10'^ Hz. For 
higher frequencies, Comptonization is thus not very efficient. 
The emergent flux is thus roughly given by the modified black- 
body spectrum corresponding to the local temperature at the ef- 
fective depth of formation. As we see in Figure 14, this temper- 
ature, which is proportional to Q, is decreased by the influence 

1 7 

of metals for ~ 2 x 10 , and so is the emergent flux. 



13 



46 



44 



42 



g> 40 



38 



10 





1 1 1 


1 


















'" r V \ 
: -.1 n V 

■ I M^ 




1 1 1 


; -M \ b 

■ I'M *■ 1 

'" 1 M M ' 
M M i r 

■ \ ■' 1 ' r 

1 ■■. 



13 



10 



14 



,15 



,16 



10^" 10* 
frequency (Hz) 



10 



17 



10 



18 



Fig. 13- Integrated spectral energy distribution vL^ for three repre- 
sentative disks for M = 10*, lO'', and 10* Mq, L/Lsad « 0.3, a = 0.01. 
The solid lines show models with metals, while the dashed Unes show 
H-He models. Dotted lines are the integrated energy distributions as- 
suming the black-body flux for each annulus. The curves are labeled 
by the values of logM. 
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Fig. 14- Computed values of AQ{T* \ x, and the Compton parame- 
ter Jeff for the annulus at R/Rg = 1.5 of the disk with M = 10* Mq, 
i/LEdd ~ 0.3. and a — 0.01. Thick lines represent the values for the 
model with metals, while the thin lines represent the values for the H- 
He model. The influence of metal edges upon the thermal destruction 
parameter e and thus yes are clearly seen. As is illustrated in Figure 16, 
the discontinuities in the curve of 4&{T*) correspond (with increasing 
frequency) to the Lyman edges of H I, He H, C VI, N VII, O VHI, 
NeX, and Fe XXIV. 

In order to show the effects of metals more clearly, we plot in 
Figure 15 the ratio of the integrated flux for models with metals 
to those for the H-He composition, for the same disks as shown 
in Figure 13. The reason for a different shape of the curves for 
different black hole masses is easy to understand. In the present 
circumstances, the metals contribute significantly to the opacity 
only through K-shell ionization of C, N, O, and Ne, and L-shell 
ionization of Fe. Only those species with edges well above kT 
retain significant abundance. Consequently, as the characteris- 
tic temperature rises, the mix of heavy elements contributing to 
the high-frequency opacity shifts to higher atomic number and 
therefore opacity at higher frequencies. All edges appear blue- 
shifted in the integrated spectrum due to relativistic effects. 
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Fig. 1 5- The ratio of integrated spectral energy distribution for the 
model with metals to that for the H-He model for three representative 
disks with M = 10*, lO', and 10* Mq, L/Lem « 0.3, and a = 0.01. 
The curves are labeled by the values of logM. We also show the posi- 
tions of the photoionization edges of the hydrogenic and He-like ions 
of He, C, N, O, and Ne, and of various ions of iron. 

The effects of changing the value of the viscosity parame- 
ter a are displayed in Figure 16. Here we display the same 
models as in Figure 13, with metals, together with analogous 
models for a = 0.1. As before, we plot the black-body energy 
distribution as well. The spectral energy distributions coincide 
for 1/ < 10'^ Hz; therefore only the high energy portions of the 
spectra are shown. 
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Fig. 16- Integrated spectral energy distributions uLi, for three repre- 
sentative disks with M = 10*, 10^, and 10* Mq, L/Ltdd ~ 0.3, com- 
puted with the influence of metal continuum opacities. Solid lines 
show models with a = 0.01, while dashed lines show models with 
a = 0.1. Dotted lines are the integrated energy distributions assuming 
the black-body flux for each annulus. The curves are labeled by the 
values of logM. 

Models with a = 0.1 generally exhibit a lower flux at fre- 
quencies below and around the maximum of flux, while the flux 
for high a models is much larger for the high-frequency region. 
This behavior is easily understood from the discussion of the 
individual annuU presented in the Sect. 4.2. Interestingly, the 
flux for a model with a = 0. 1 for a given M is very close to the 
flux predicted for a model with a = 0.01 and M about an order 
of magnitude lower! We also conclude that for the hot disks, at 
least those considered in the present paper, the effects of chang- 
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ing the value of the viscosity parameter a are more important 
than the effects of metal continuum opacities. 

Finally, as a reference point for later work that examines 
surface ionization response to different heat deposition laws, 
we show the ionization structure for three representative an- 
nuli of the disk with M — 10'' M0, namely the hottest annu- 
lus, R/Rg = 1.5 with Teff = 807,000 K; an intermediate annu- 
lus, R/Rg = 15 with Jeff = 207,000 K; and a "cool" aimulus, 
R/Rg = 150 with Teff = 40,000 K. We show the ionization for 
two species: oxygen, which is a representative for the light met- 
als; and for iron. The knowledge of ionization balance is im- 
portant to show us what are the most important opacity sources 
in the X-ray region. The results are shown in Figure 17 for 
oxygen, and Figure 18 for iron. The curves are labeled by the 
effective charge of the ion, i.e., "1" denoting neutrals, "2" once 
ionized, etc. The label appears close to the position of the max- 
imum number density of the given ion. 

Oxygen is essentially fully stripped in the hottest annulus 
throughout the whole vertical extent; while He-like oxygen (O 
VII) dominates around r w 1 for the intermediate annulus. In 
the cool annulus, O 111 and IV are dominant stages around a 
Rosseland mean optical depth r 1, while O V and VI domi- 
nate deeper down. The ionization of iron is shown in Figure 18. 
The dominant stages of iron are Fe XXll - XXV for the hot an- 
nulus; Fe XV - XVII for the intermediate annulus, and Fe IV - 
VII for the cool annulus. 
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Fig. 17- Ionization structure of oxygen for three representative annuli 
of a disk model with M = lO** M0, L/LEdd « 0.3, and a = 0.01, at 
radial distances R/Rg = 1.5, 15, and 150. These annuli have effective 
temperatures approximately 807,000; 207,000; and 40,000 K, respec- 
tively. The abscissa is the column mass in g cm^" and the ordinate 
is the ion number density in cm~^. The curves are labeled by the ef- 
fective charge of the ion, i.e., "1" denoting neutrals, "2" once ionized, 
etc. 
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Fig. 1 8- The same as in Figure 17, but for iron. 

5. DISCUSSION 
5.1. Effect of Comptonization on H and He Edges 

The observed absence of edge features at the hydrogen Ly- 
man limit has been a commonly cited problem with accretion 
disk models of AGN. As we discussed extensively in Paper III, 
our disk models do not exhibit such features for accretion rates 
which are reasonable fractions of the Eddington Umit, at least 
for the high black hole masses considered there. Comptoniza- 
tion only reinforces this conclusion, because it tends to smear 
out flux discontinuities (Czerny & Zbyszewska 1991). This is 
clear in the spectra of individual annuli shown in Figure 7. The 
effect is most dramatic for the Hell Lyman edge, which occurs 
at the high frequencies where Comptonization is most impor- 
tant, making it a much smoother feature (cf. Figure 8). 

Our integrated disk spectra in Figure 1 3 show that the H and 
He II Lyman edges become more pronounced as the black hole 
mass is lowered. This effect is easily understood. The vertically 
averaged density as a function of R/Rg scales as (see Paper II) 
poc {R/Rgfl^O{R/Rg)MM-^a-\ where the factor 0(R/Rg) 
accounts for relativistic corrections. Since we consider mod- 
els with fixed L/L^ii, we have L/LEdd oc M/M = const, i.e., 
M cx M, so that the density scales for our set of models as 
p (X M~^a~^. The models with lower mass have higher den- 
sity, and therefore higher thermal coupling parameter e, so that 
the bound-free edges are stronger. The above considerations 
also indicate that the edges are weaker for higher values of a 
(cf. Figure 16). Our results suggest that it might be easier to 
find broadened edge features in Seyferts rather than quasars, al- 
though one would have to look out to reasonable redshifts in 
order to avoid Galactic interstellar absorption. 

5.2. Metal Opacities and Soft X-ray Emission 

In these disks, whose heating rate per unit mass is either con- 
stant or declines toward the disk surface, the prime effect of 
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heavy element K-shell opacity is to partially suppress emission 
in the soft X-ray bands affected by these photoionization edges. 
The particular edges that most powerfully affect the spectrum 
depend on the central mass, rising in energy as the central mass 
decreases (Figure 15). To a very rough approximation, the cen- 
tral energy of the edge that produces the greatest X-ray absorp- 
tion scales as M"'/"*, similarly to the temperature scale in the 
disk. As a result of this additional opacity, flux is redirected to 
frequency bands with smaller opacity in the EUV (Figure 13). 

However, there is reason to think this effect may be strongly 
model-dependent. As shown in Figure 16, its magnitude can be 
substantially altered by changes in the local column density (as 
parameterized, for example, by the stress parameter a). Thin- 
ner disks present less soft X-ray opacity, and therefore show 
weaker breaks in that region than do thicker disks. Similarly, if 
the heating rate per unit mass were to increase, rather than de- 
crease, in the surface layers, these same features might swing 
into emission. 

Observations of quasar spectra in the ^ 100 — 1000 eV 
band reveal a complicated (and possibly contradictory) story. 
Low redshift radio-quiet quasars appear to be consistent with 
single steep power-laws with an average spectral slope of 
<ilnL^/<ilni/ ~ -1.7 (Laor et al. 1997). On the other hand, 
in many cases the spectra of lower luminosity nearby radio 
quiet AGN, while following power laws above ~ 2 keV, rise 
sharply at lower energies (Turner & Pounds 1989), producing a 
so-called "soft excess". Similarly, Blair et al. (2000) find that 
the composite X-ray spectra of quasars binned over various red- 
shift intervals are best fit by a combination of a power-law and 
single-temperature blackbody. The power law component has 
a spectral slope of d\nLi, / d\nv ~ —0.9, consistent with the 2- 
10 keV slope seen in nearby Seyferts. The blackbody compo- 
nent that they fit has a temperature that rises from ~ 100 eV at 
z ~ to ~ 270 eV at z ~ 2.5. The luminosity of this component 
is ~ 10"*"^ ergs s~', independent of redshift. 

As shown in Figure 16, the spectral energy distributions ex- 
pected from our disk models in the 0.2-2 keV (5 x lO'^ — 5 x 
10'^ Hz) band show a broad range of slopes. At the highest en- 
ergies they all show exponential rollovers, however, so the hard 
power-law spectra observed in AGN require a separate compo- 
nent. This is widely believed to be due to thermal Comptoniza- 
tion in a hot corona. A modification of our models to include 
increased viscous dissipation near the surface might be able to 
self-consistently produce such a corona, and we will explore 
this in a future paper 

Whether or not the models presented in this paper can pro- 
vide good fits to the soft excess alone is unclear Our models 
require L/LEdd < 0.3 in order to be consistent with our assump- 
tion that the disk is geometrically thin. A black hole mass of at 
least 10^ M0 is required to produce observed quasar luminosi- 
ties, and Figure 16 shows that our accretion disk models for 
such black holes will produce spectra that peak below 200 eV. 
This might nevertheless be consistent with the 270 eV tempera- 
tures of the Blair et al. blackbody components at high redshift. 
Our models are not well-fit by single temperature blackbodies 
near the peak, so detailed spectral fits are required. 

Huge soft X-ray excesses are a hallmark of narrow line 
Seyfert 1 galaxies (Boiler, Brandt, & Fink 1996). In order to ex- 
plain the soft X-rays, accretion disk models of these sources are 
inevitably forced to high accretion luminosities (Eddington or 
above), if they are based on the assumption of local blackbody 
emission (e.g. Mineshige et al. 2000). However, our models 
shown in Figure 16 indicate that these high accretion rates may 



not in fact be necessary. Comptonization extends the emission 
well beyond that expected from a multi-temperature blackbody 
model of the same disk. This might contribute substantially to 
the observed soft X-ray excesses in narrow line Seyfert Is, par- 
ticularly at the relatively low black hole masses that might exist 
in these sources. 

5.3. Radiative vs. Convective Heat Transport 

We have assumed throughout the present paper that heat is 
transported vertically through the disk by radiative transfer. 
However, it turns out that all our disk models are convectively 
unstable. Away from ionization zones, the square of the Brunt- 
Vaisala frequency in the optically thick limit may be written 
as 



dlnp 



3P 



dz V5^'gas+4Prad 

where the adiabatic gradient is 
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and P = f,ad + ^'gas is the total (radiation plus gas) pressure. 
Note that the adiabatic gradient is 4/3 for a radiation pres- 
sure dominated medium, but 5/3 for a gas pressure dominated 
medium. Equation (^) shows that convective instability oc- 
curs when dlnP/dlnp exceeds the adiabatic gradient, with a 
characteristic growth rate of order the orbital frequency. Fig- 
ure 19 depicts lilnP/iilnp as a function of column density for 
a number of annuli in ourM = 10^ M0, L = 0.3 Lem, ol = 0.01 
model. We see that all our annuli are convectively unstable at 
high column densities. Photon diffusion is not included in equa- 
tion (^), but the effect of a long photon mean free path near the 
photosphere would be to diminish the radiation pressure con- 
tribution to the adiabatic gradient, which therefore approaches 
5/3. 
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Fig. 19- Plot of the logarithmic gradient dInP/dlnp as a func- 
tion of column density for a number of annuli in the M = 10* M0, 
L = 0.3 LEdd, a — 0.01 model. The curves are labeled by the value of 
R/Rg- Dotted lines show the adiabatic gradient for the hottest annulus 
(R/Rg = 1.5 - lower curve), and for the coolest annulus (R/Rg = 100 
- upper curve). In both cases, the adiabatic gradient is close to 4/3, al- 
though it is higher for the cooler model, where the gas pressure makes 
a non-negligible contribution. 

These instabilities depend of course on our assumed verti- 
cal dissipation profile. At the present time this profile is un- 
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known, although numerical simulations of the nonlinear devel- 
opment of the magneto-rotational instabihty are beginning to 
address this question, at least in gas pressure dominated disks 
(e.g. Miller & Stone 2000). It is nevertheless clear that our 
models are not self-consistent, and convective heat transport 
should be included. How this will affect the spectrum is unclear 
While the scattering photospheres appear to be convectively 
stable, the continuum is formed in deeper, potentially unsta- 
ble layers. More seriously, hydrodynamical simulations show 
that convective heat transport changes the overall heat content 
in the disk, altering the scale height and surface gravity of the 
photosphere (Agolet al. 2000). These simulations also produce 
an approximately isentropic configuration (at least when hori- 
zontally averaged), similar to what is seen in convection zones 
in stars. Assuming that magnetic fields do not change things 
(see, however, Hawley, Balbus, & Stone 2001), the assumption 
of constant entropy may provide a good first approximation to 
convective heat transport in unstable zones. We will address 
convective heat transport in disks in a future paper. 

6. CONCLUSIONS 

We have extended our modeling procedure described in the 
previous papers of this series in two respects. First, we have 
self-consistently included the effects of Compton scattering, 
both on the disk structure and on the emergent radiation. Sec- 
ond, we have included the effects of continuum opacities of the 
most important metals (C, N, O, Ne, Mg, Si, S, Ar, Ca, Fe, Ni). 

The present models are thus based on the following assump- 
tions: 

(i) The energy is generated by turbulent viscous dissipation, 
with the vertically-averaged viscosity described through the 
Shakura-Sunyaev parameter a. 

(ii) The kinematic viscosity is constant with depth. 

(iii) Convection and conduction are neglected. 

(iv) No external irradiation or self-irradiation of the disk is con- 
sidered. 

(v) Effects of line opacity are neglected. 

The underlying assumption, not listed here, is that the 1-D ap- 
proach is appropriate, i.e., that the disk may be described as a 
set of mutually non-interacting, concentric annuli. 

We have also developed a simple analytic model that pro- 
vides an estimate of the vertical temperature structure of the 
disk in the presence of Comptonization. An interesting by- 
product of this model is an estimate of the surface temperature. 
In particular, we have shown under which conditions a hot outer 
layer (a "corona") is formed. 

We have recomputed a number of models from the very large 
central-mass grid of Hubeny et al. 2000 (Paper III) with Comp- 
tonization taken into account, and showed that the predicted 
disk-integrated spectra differ only modestly from the previous 



models. We have concluded that the grid of models from Pa- 
per 111 is not significantly plagued by neglecting Comptoniza- 
tion, except perhaps for the very hottest (smallest central mass) 
models of the sample. We have also compared our results to 
analogous models of Ross, Fabian, & Mineshige (1992). The 
agreement of the computed continuum flux is very good, al- 
though there are some differences, particularly at the shortest 
and the longest wavelengths. The predicted profiles for the hy- 
drogen and helium lines differ significantly, which is very likely 
explained by differences in treating the lines, in particular by 
the escape probability treatment adopted by Ross et al. 

We have also computed a set of vertical structure models for 
black hole masses of 10'', 10^, and 10** M0, total luminosity 
L/Lem ~ 0.3, and for two values of the a parameter: 0.1, and 
0.01. Unlike the previous papers of this series, we consider 
a depth-independent kinematic viscosity. The main effect of 
Comptonization in this mass range is to transfer about 10 per- 
cent of the total luminosity from the high energy tail (above 
6 X 10'^ Hz) to the peak of the spectrum (between 1.5 and 
4 X 10'^ Hz.) We have computed a number of tests, and found 
that for the physical conditions studied in this paper (tempera- 
tures up to 10^ K), the effects of Mg, Si, Ca, and Ni are negli- 
gible. We have shown that the effects of other metals are still 
relatively modest. The only significant differences in the pre- 
dicted spectra occur in the high-frequency tails of the spectra, 
where the flux is already rather low. Interestingly, the effects of 
changing the viscosity parameter a from 0.01 to 0.1 are found 
to be larger that the effects of metals. Models with a = 0.1 
generally exhibit a lower flux at frequencies below and around 
the maximum of flux, while the flux for high a models is much 
larger for the high-frequency region. 

The effects studied in this paper are most dramatic in the 
EUV/soft X-rays. As shown in Figure 16, this high frequency 
region of the spectrum differs enormously from that expected 
from multitemperature blackbodies. Fits to the soft X-ray ex- 
cesses of AGN which assume local blackbodies will necessarily 
require an artificially high accretion rate and accretion luminos- 
ity. 

In future papers of this series, we will relax assumption (ii), 
and will compute models that allow for an increase of the heat- 
ing rate (kinematic viscosity) with height, leading to the for- 
mation of very hot upper layers of the disk. We will also relax 
assumptions (iii), incorporating convective heat transport, and 
(v), studying the effects of line opacities on the disk structure 
in detail. 

We thank Randy Ross for providing us with his model in dig- 
ital form. We also thank Tim Kallman for sharing with us his 
collection of atomic data. This work was supported in part by 
NASA grant NAG5-7075. 



APPENDIX 

APPENDIX A: THE COMPTON SCATTERING SOURCE FUNCTION 

We treat Comptonization in the non-relativistic, diffusion limit, where the exact integral over the redistribution function is approx- 
imated by a Kompaneets type operator. The general form of the radiative transfer equation for arbitrary geometries in this limit has 
been derived by Babuel-Peyrissac & Rouvillois (1969) and Nagirner & Poutanen (1993).[| In the plane-parallel geometry considered 
here, the radiation field is axisymmetric, but still has a dependence on the polar direction cosine After integrating over azimuthal 
angles, we find that the contribution of Compton scattering to the right-hand side of the radiative transfer equation (i.e. sources minus 

^The equations in these papers contain small typos which we have connected. 
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sinks) is 



Compt 



^n,CTT J dm [{Ii-I){A-2B + C)+B{D+E)] . 



(Al) 



For simplicity of notation here, we have dropped the frequency subscript, u. The specific intensity / = while h = IvifJ-i)- 
Here, 

A(/i,/ii) = a, (A2) 



B{n,iJ,i) =x{a-l3), 
C{n, Ml )= 26 (8 - 87 - 3a + 2/3) , 



D = 



and 



where 



E = A 



h dv Xv^ 



-h 1 



2hv^ 



hv 



X = 



6 = 



kT 



(A3) 
(A4) 

(A5) 



(A6) 



(A7) 



nteC^ meC- 
and where the matrices a, (i, and 7 are given by 

a{^J.,^Jil) = ^^?^i\-^?-^l\+^, (A8) 

0{lJ,,lii)=5iJ,iJ,i+5/?/4-3/?/ii-3nfj,l, (A9) 
7(M,Mi)=m, (AlO) 

so that a is an even function of fi 01 fii, while /3 and 7 are odd. 

The radiative transfer equation is conveniently written by splitting the total absorption/emission term into four parts: absorption 
{(^vlv), spontaneous enussion {r}v), stimulated emission {rf^^Iv), and thermal absorption/emission (k*, 77*), i.e.. 



where 



3 

a^ = —ne(7Tj diJbi{A-2B + C), 
'n^ = ^ne<Ji j dfii[h{A-2B + C)+B{D' + E% 



3 r 
_stim _ 



2 ri 



dniB 



2hv^kT d f h 
kT h dv \ 



-4/1 



where D' and £" are the D and £ terms without the stimulated enussion contributions, i.e.. 



D' 



kT d /I 



yp- dv^ j h dv \ 



and 



We introduce the following notation 



so the D' and E' terms are written 



£' = 4 



kT ,9 ( h 



h dv Kv^ 



I' 



I" 



+ h 



and thus 



^ dlnv' " dilnvy 



kT 

D' = ^ {!{' - 1I[ + 12/1) + (/( - 3/1) , 

kT 

E'=A-{l[-3h)+AIu 



x{D' + E') = QI[' + (x- 39)/; +xli . 



(All) 

(A12) 
(A13) 
(A14) 

(A15) 
(A16) 

(A17) 

(A18) 
(A19) 

(A20) 
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After some algebra, we obtain for the individual terms 

= Heaj (1 — 2x) , 



(A21) 



and 



r]^ = ^neaT [ <iAii{ [(1 -x- 6e)a+ 166 + (x + 4e)/3- 1667] 7i 

+ {a-(3){x-3e)l[ + {a-f3)Ql['}, (All) 

??r = ^«.^T^ j dni2x{a-P){l[~h). (A23) 

The thermal terms k* and ry* contain bound-free, free-free, and possibly also bound-bound contributions from all transitions which 
are taken into account. The stimulated emission is treated as negative absorption, as is customary in astrophysical radiative transfer. 

The Compton emission terms are more complicated than the thermal emission term r/* because they depend explicitly on angle. 
To write them in a more useful form, let us first introduce the traditional moments of the specific intensity, defined by 

[J,H,K,N]=- J ^dnih-[l,nufA,fA]- (A24) 
and also special "pseudo-moments" of the specific intensity with weights a, (3, and 7, defined as 

VJa, Jp. Jj]=^J dinh ■ [a, /?, 7], (A25) 
The pseudo-moments can be expressed through the ordinary moments as 

Ja = \[{'i-IJ^)J+{2>IJ^-l)K\, (A26) 

^/3 = ^m[(5-3/x2)H+(5m'-3)A^], (A27) 

J-, = \lJiH. (A28) 

Because the pseudo-moments are integrals over matrices a, (3, and 7, they are exphcit functions of the angle, /i. We may then also 
introduce moments of Ja, J(3, and J-^. The zero-order moments are given by 



,dfi = J, (A29) 

and (since jT/j and are antisymmetric functions of fj,), 

jpdiJ.= -l J^diJL = 0. (A30) 



1 . 1 



The first-order moments are given by 



and 



1 

JaHdiJ, = 0, (A31) 

■ 1 



JpHdn='^H, J^fidn=^H. (A32) 

Using the definitions of the optical depth, t^, the photon destruction probabihty, e^, and the scattering probabihty, X^, introduced 
in section 2.1, we obtain for the full, angle-dependent transfer equation 

^dl^ = - eX - A.5^°'"P'(/^) , (A33) 
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where the (generally angle-dependent) "Compton source function" is given by 

+Ufi)b,2x[j;^{fi)-Jf,{fi)-J^{ti) + Jis{fi)], (A34) 
where the last term represents a contribution of the stimulated Compton scattering, with 

^2 

(A35) 



2/11/3 



So far we have not made any approximatio ns bey ond the Kompaneets limit of low photon energies (jc <C 1) and low electron 
temperatures (0 ^ 1). In particular, equation (\34) retains the full angle dependence of the Compton scattering source function. 
However, it is quite complicated to solve numerically, and we have made further approximations for all the calculations presented in 
this paper. 

The first approximation deals with the stimulated emission term, which is nonlinear in the specific intensity. Fortunately, it is 
usually very small compared to the spontaneous term in our models, and could be neglected completely. We choose to retain it, and 
show how to treat the nonlinearity below. However, we approximate the angular dependence of the radiation field by replacing 
with Ji, in the stimulated term. 

More g enera lly, we average the angle de pend ence of the Compton scattering source function by first taking moments of the transfer 
equation ( ^33 ). Using equations ( A.26) - (A32), the first and second moments are 

r)H 

—^=J,- e,5* -K[{l-x)J, + {x- 3e)Jl + ef^] - KJ,b,,2x{Jl - J,) , (A36) 
dK 2 2 

-^=H,-K-[{x-Q)H,-{x- 3Q)Hl - QH'J] - A. - J,bax{H, - Hi) . (A37) 

We now assume that all terms on the right-hand side of equation ( A37 ) are negligible with respect to the leading term H^, so that the 
second moment equation is 

^=H.- (A38) 

We now introduce the variable Eddington factor, 

U=K,/J,, (A39) 
so that the two moment equations are combined to a single equation. 



drl 



J, - e,5j - A,[(l - x)J, -^{x- 36)4 + " " ■ (A40) 



This equation is analogous to equations (|I]) and (^, and is used for computing the disk structure. 

We introduce a set of discretized frequency points, i'i,i=l,...,N; with ly, < Vjj^ i . We denote /,■ =I{vi). The frequency derivatives 
are treated using the following discrete representation. Let us first take internal frequency points, i = 2, . . . ,N — I. Then 

//=cr/,._i+c^/,. + c+/,+i, (A41) 

//' - Jr/,_ 1 + d,"/,. + d+I^+i , (A42) 
where the coefficients , cf, c'^ are determined by the method of Chang & Cooper (1970), and 

^r = ^^, (A43) 

d^ = ^^' (A44) 

d^ = -dr^d+, (A45) 

where 

A/_i/2 = ln(z/,-/i/,-_i) , A,+ i/2 = ln(j/,-+i/t/,-) , (A46) 

and 

A,- = A;_i/2 + A/+i/2. (A47) 
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The boundary conditions in frequency space are written in the form of equations ( A,41 ) and ( A42 ), with Cj" = — — — 0, 
and c'^ = d^ = c'^ = d^ = 0. At the lowest frequency we have, assuming the Rayleigh-Jeans form of the specific intensity, equation 

c0 = 2, d\=A. (A48) 

At the highest frequency we have either 

c°=3-z, < = (3-z)2-z, (A49) 
(with z = hu/kT) if we assume the Wien form of the specific intensity, equation (|l^; or 

4 = -/3, dl = l3\ (A50) 

if we assume the power-law form, equation (p^. 

We write the radiative transfer equation for the /-th frequency point as 



drf 

where 



= Ji - eiSf - A,[ (A,y,_ 1 + BiJi + C,J,+i ) - /,■ (£,/,■ + i + y,J,+ i) ] , (A5 1) 



Bi = (1 -Xi) + [xi - 3e)c[' + e^/f , (A52) 

Ai = {xi - 3e)cr + Qdy , (A53) 

Q = [xi - 3 e )c+ + Qd+ , ( A54) 

£,=/7/2x,(c°-l), (A55) 

Ui = bilxicr , (A56) 

Vi = bi2xic+ . (A57) 

We now avoid the nonlinearity of the stimulated emission term by evaluating it partly by using "old" mean intensities. In other words, 
we replace B, by where 

b; = Bi + Eijf^ + UiJf}, + ViJf^, . (A58) 

The discretized transfer equation is then written 

= Ji - eiSf - Xi {AiJi_ 1 + B'iJi + C,y,+i ) , (A59) 
The final step is to consider a discrete representation of the differential term as in the standard Feautrier method (Mihalas 1978), viz. 



dr} 



aijJi4- 1 - PijJi^ + jijJi,d+ 1 , ( A60) 



where index d labels depth points, while / labels the frequency points. The coefficients a, (3, and 7 (not to be confused with the 
previously introduced angle-dependent matrices of Compton scattering) are given by 

ai4 = //,rf-i/(Ar,-.j_i/2Ar,y) , (A61) 

li.d = fi.d+\l (At,- rf+i/2Ar;.d) , (A62) 

= Il!L ( ^ + ^ ^ , (A63) 



Ar/,rf \Ar,'rf_i/2 Ar,-,rf+i/2 
and 

^Tij_ I /2 = Tij - Tij^ 1 , ( A64) 

Ar/.rf+i/2 = Tij+i - Tij , (A65) 

Ar;,rf = (t/,j+i/2 + T-,-,d-i/2)/2. (A66) 
The (second-order) boundary conditions are given by 

Am = 2/;« /At,3/2 + 2/;- 1 /( At,3/2)2 , (A67) 

7,-,i =2/,■,2/(AT,3/2)^ (A68) 
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and 



A-,D = 2Az,/(At,b_,/2)2, 

ai,D = 2//.D_i/(Ar,-£,-i/2)^ • 
The final discretized form of the transfer equation thus reads 

— aijJij-i — ji,dJi,d+\ + {Pi,d + 1 — XijB'jj)Jij — XijAijJi-ij — \i.dCi,dJi+\,d = £i.dS% ■ 



(A69) 
(A70) 

(A71) 



The formal solution, i.e., a solution of the transfer equation with the known thermal source function, 5*, is obtained as follows. We 
organize the mean intensities in a set of column vectors 



J/ = (Ji,l,Ji,2, ■ ■ ■ ,Ji,Dj 



(All) 



so that the mean intensity vector contains intensities in all depth points for a given frequency point /. The resulting discretized transfer 
equation reads 

- A,J,._ 1 + B,J,- - C,J/+i = U , (A73) 

where the matrices B are tridiagonal (because of the difference representations of the second derivative with respect to depth), while 
the matrices A and C are diagonal (because the terms containing the frequency derivatives are local in physical space). 
The matrix elements are given by 

(A74) 



{^i)d.d = Pi.d + 1 — '^i.dB'i d I 
(^i)d.d-\ = —Ctij , 
{^i)d.d+] = —li.d , 
{^i)d.d = \.dAij , 
iCi)d.d = Xi.dCij , 

(L,)j = <^i,dSf^d- 



(A75) 
(A76) 
(A77) 
(A78) 
(A79) 



Solution of equation ( 473 ) is done by the standard Gauss-Jordan elimination, consisting of a forward elimination followed by a 
backsubstitution. 



D, = (B/-AiD,- 



i = 2, 



with 

and 

with 

The back substitution step is 



Di=B-'Ci, 



Z, = (B,--AA-i)"' (L,- + A,-Z,-_,), 
Zi 



:B-'Ll 



J,=D,J, 



+ 1 ■ 



(A80) 
(A81) 
(A82) 
(A83) 
(A84) 



with Jn = Zn- Starting at / = 1, we compute successive values for and Z, through i = N — 1. At the last point, / =JV, Cn = 
and hence Dai = 0, and ~ Z^. Having found Jn, we then perform successive back-substitutions into equation (/V84) to find 
J,,/ = A^-1,...,1. 

APPENDIX B: TREATMENT OF DIELECTRONIC RECOMBINATION 
In terms of a cross-section, the photoionization rate is given by (see, e.g., Mihalas 1978) 



f°°4TT , ^ 

Rik = / —(j[v)J^dv , 



(Bl) 



where <t{v) is the photoionization cross-section, v the frequency, v^) the frequency of the continuum edge, and J,j the mean intensity 
of radiation. The subscript ; labels the starting level of the bound-free transition, while k labels the ending level (usually the ground 
state of the next higher ion). 

The radiative recombination rate is given through detailed balance arguments (the Einstein-Milne relations for the continuum; see 
also Mihalas 1978) as 



f°° 47r flhv^ \ 
Rki = n^(j)i{T) J —a{v)y—^+J^jtx^{~hv/kT)dv, 



(B2) 
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where «e is the electron density, and (j>i{T) the Saha-Boltzmann factor, given by (Mihalas 1978, p. 113) 

/ ,2 \ 3/2 

<t>i{T) = ] ^ cxp{h,,o/kT) = 8i exp(/,j.o/^r) , (B3) 

yimnkT J 2gk gk 



where g are the statistical weights. The first term on the right-hand side of equation (B2) represents spontaneous recombination, 
while the second term (proportional to J^) represents stimulated recombination. 

We now use this to express the known dielectronic recombination rate through an artificial cross-section. For simplicity, we assume 
that 

= (To , for < v < y\ , (B4) 

and is zero everywhere else. We further assume that vi/v(_) — \ ^ 1, i.e., the frequency v\ is not very different from vq. We further 
assume that the stimulated term is negligible (i.e., J^, ^ 2hv^ /c^). The known recombination rate is given through the artificial 
cross-section as follows ^ 

R = n^CiT^^^^ — exp{hi^o/kT) [ —ao'^^exp{-hv/kT), (B5) 
gk Juo hv 

or 

R = Cj^n^T-^/^ — ao / ly^ exp[~h{i^ - i^o)/kT] . (B6) 



'1^0 

gk 

The integrand is the product of two terms. The first one, v^, varies slowly with v compared to the second, exponential, term. We may 
therefore remove the term from the integral, replacing it with = {{vq + v\) /2]^ ~ Vq, and performing the remaining integration 
analytically, 

R = C!^n^T-y^^aov^\\-exp{-~hl:\vlkT)\ , (B7) 
hc^ gk 

where /S.v — v^ — vq. In the calculations reported in this paper, we used = 1 . 1 1'o- 
Consequently, the desired value of the cross-section, co, is given by 



ao = R— Q ^ - 4 [1 - exp {-h/^^v/kT)] ' . (B8) 
He oi:k gi v 



Unlike the traditional cross-section, cro now generally depends on temperature and electron density. We treat this dependence exactly. 
Nevertheless, the dependence is rather weak because the dielectronic recombination rate is typically proportional to T^'/^, so the 
T-dependence is contained essentially in the term exp[— — v^j/kT}, which is typically small compared to unity. Analogously, 
the dielectronic recombination rate is directly proportional to Wg, so the dependence of cto on physical parameters is to a large extent 
factored out, and one is left with ctq essentially independent of state parameters. 
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